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ABSTRACT 

Formation of bodies near the deuterium-burning limit is considered by detailed numerical simulations 
according to the core-nucleated giant planet accretion scenario. The objects, with heavy-element cores 
in the range 5-30 M$, arc assumed to accrete gas up to final masses of 10-15 Jupiter masses (Mj up ). 
After the formation process, which lasts 1-5 Myr and which ends with a 'cold-start', low-entropy 
configuration, the bodies evolve at constant mass up to an age of several Gyr. Deuterium burning 
via proton capture is included in the calculation, and we determined the mass, M50, above which 
more than 50% of the initial deuterium is burned. This often-quoted borderline between giant planets 
and brown dwarfs is found to depend only slightly on parameters, such as core mass, stellar mass, 
formation location, solid surface density in the protoplanetary disk, disk viscosity, and dust opacity. 
The values for M50 fall in the range 11.6-13.6 Mj up , in agreement with previous determinations that 
do not take the formation process into account. For a given opacity law during the formation process, 
objects with higher core masses form more quickly. The result is higher entropy in the envelope 
at the completion of accretion, yielding lower values of M50. For masses above M50, during the 
deuterium-burning phase, objects expand and increase in luminosity by 1 to 3 orders of magnitude. 
Evolutionary tracks in the luminosity-versus-time diagram are compared with the observed position 
of the companion to Beta Pictoris. 
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1. INTRODUCTION 

There is considerable debate over the question of defin- 
ing a precise boundary between the class of objects 
known as 'planets' and those known as 'brown dwarfs'. 
It has been suggested that the two types of objects could 
be distinguished by their formation mechanism; however, 
it is generally difficult to deduce this property from ob- 
servations of specific objects. Nevertheless, there is a 
well-defined minimum in the mass distribution (actually 
Msini), for substellar companions to G and K main- 
sequence stars, in the range 20-30 Mj up (Lovis ct al. 
2006; Sahlmann et al. 2011; Schneider et al. 2011), sug- 
gesting that objects near the deuterium-burning limit 
can 'form like planets'. These authors suggest that this 
minimum does correspond to a (somewhat imprecise) 
dividing line between formation mechanisms and that 
the upper limit to planet masses should be set at about 
Msini rj 25 Mj up . However, no break is seen near this 
mass in the distribution of free-floating objects observed 
in the Sigma Orionis young cluster (Pena Ramirez et al. 
2012) down to 4 Mj up . Together, the observations imply 
that formation mechanisms do not define a unique mass 
boundary between planets and brown dwarfs. 
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Another commonly used criterion to classify planets, 
brown dwarfs and stars is based on nuclear fusion that 
does or does not occur within the object. Brown dwarfs 
are defined to be those objects that at some point in their 
evolution become hot enough in their interiors to burn 
a majority of the deuterium that was initially present 
in the object; however, they never become hot enough 
to burn -"-H by the proton-plus- proton reaction in a self- 
sustaining manner as true stars do. On the other hand, 
the term planet is applied only to objects that will not 
burn much deuterium. This criterion was used by Bur- 
rows et al. (1997) to separate the two types of objects, 
and the dividing line was stated to be ~ 13 Mj up , where 
M Jup = 1.898 x 10 30 g. This dividing line depends on 
the helium mass fraction, the deuterium abundance, and 
the metallicity, and Spiegel et al. (2011) found that for 
a reasonable range of parameters, 50% of the initial D 
is burned in the mass range 12-14 Mj up . The evolution- 
ary models used to establish the criterion have a uniform 
chemical composition, a defined total mass in the vicinity 
of the 13 Mj up limit, constant in time, an initial radius 
of about 2-3 Rj up where Rj up = 7.15 x 10 9 cm, and an 
initial photosphcric temperature (T e g) of about 2500 K. 
The corresponding initial luminosities are 2— 3 x 10~ 3 L Q . 
A starting model of this type has become known as a 
'hot-start' model, characterized by a relatively high ini- 
tial entropy (Marley et al. 2007). 

The question of whether, either for brown dwarfs 
or planets, the formation mechanism actually leads to 
such hot-start initial conditions is still under investiga- 
tion. For objects formed either by collapse of interstellar 
clouds or by fragmentation in a protostcllar disk by grav- 
itational instability, it is plausible that the hot-start ini- 
tial condition could be reached (Baraffc et al. 2002). In 
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the case of gravitational instability, Galvagni et al. (2012) 
have found, from three-dimensional numerical simula- 
tions, that the entropy of newly-formed clumps, near the 
point where molecular dissociation sets in at the center, 
is high, possibly consistent with a hot start. However 
the dominant process for giant planet formation is most 
likely the core-nucleated accretion mechanism, in which 
solid particles first accumulate to form a heavy-element 
core, then later when the core has attained roughly 5- 
10 M®, gas is captured from the disk. A particular set of 
evolutionary calculations based on this theory (Marley 
ct al. 2007) shows that once the planet has become fully 
formed, its entropy is relatively low, with luminosities on 
the order of 10 -5 — 10~ 6 L Q . The low entropy is a direct 
consequence of the assumption made in these calcula- 
tions that, during the phase of rapid gas accretion, all 
of the accretion energy is radiated away at the accretion 
shock at the planet's surface. Thus, the core accretion 
process can lead to a 'cold start'. However the shock 
treatment is approximate, and the accretion flow cannot 
actually be modelled correctly with 1-D spherically sym- 
metric calculations. Thus other possibilities can arise. 
Mordasini et al. (2012) show that core accretion forma- 
tion calculations in which none of the energy is radiated 
at the shock lead to hot-start conditions very similar to 
those assumed by Baraffe et al. (2003) and Burrows et al. 
(1997). Furthermore, intermediate 'warm' states are also 
possible outcomes (Spiegel & Burrows 2012). In the core- 
accretion picture, also, the chemical composition is not 
uniform because of the presence of the core, which turns 
out for the case of a Jupiter mass planet to fall in the 
range 4-20 M e (Movshovitz et al. 2010). 

A massive object of 25 Mj up formed by core accre- 
tion (Baraffe et al. 2008) has been shown to burn all 
of its initial deuterium despite the presence of a heavy- 
element core of 100 or a few hundred M®. Cold-start 
models, including the core and calculations of the for- 
mation phase, have been investigated to determine the 
D-burning mass limit (Molliere & Mordasini 2012). The 
results show that the limit still falls within the range 12- 
14 Mj up . The purpose of the present paper is to present 
further formation calculations for bodies formed by core- 
nucleated accretion that end up with a total mass in the 
10-15 Mj up range in the low-entropy state, and to inves- 
tigate the effect of various possible initial conditions, as 
well as physical parameters during the formation stage, 
on the corresponding deuterium-burning limit. 



2. COMPUTATIONAL METHOD 

The evolutionary calculations for giant planets are 
started at the point where the heavy-element core has 
a mass of about 1 M®, and are carried through the en- 
tire formation process as well as the subsequent contrac- 
tion/cooling phase at constant mass, up to a final age of 
several Gyr. The assumptions and computational pro- 
cedures were described in detail in previous publications 
(Pollack et al. 1996; Bodenheimer et al. 2000; Hubickyj 
ct al. 2005; Lissaucr ct al. 2009; Movshovitz et al. 2010). 
The early phase of the formation process is dominated 
by the accretion of planetesimals onto the core; during 
this phase the gaseous envelope has low mass, <C 1 Mg, 
and a low accretion rate compared to that of the core. 



The latter is given by 

= nR 2 cllpt an p F g (1) 

where 7ri? 2 apt is the effective geometrical capture cross 
section, a is the surface density of solid particles (plan- 
etesimals) in the protoplanetary disk, Q p is the planet's 
orbital frequency, and F g is the gravitational enhance- 
ment factor, which is obtained from the calculations of 
Greenzweig & Lissauer (1992). The planetesimal radius 
is taken to be 50 km for the cases with a central star of 
2 Mq and 100 km for the cases with a star of 1 Mq (see 
Tabic 1). The smaller size, or a reasonable distribution 
of planetesimal sizes, tends to reduce the formation time 
but has little effect on the basic results of this paper. 

If no gaseous envelope is present, then i? C apt = Rcorc, 
the radius of the heavy-element core. However, even if 
the envelope mass is relatively small compared with the 
core mass, the planetesimals interact with the envelope 
gas, are slowed down by gas drag, and are subject to 
ablation and fragmentation. The trajectories of plan- 
etesimals through the envelope are calculated (Podolak 
et al. 1988), and the effective -R cap t is determined. The 
material that is deposited in the envelope is then allowed 
to sink to the core, as discussed by Pollack et al. (1996). 
Calculations by Iaroslavitz & Podolak (2007) show that 
this assumption is valid at least for the organic and rock 
components of the planetesimals. The ices, however, can 
dissolve in the envelope, so that our 'core mass' is some- 
what overestimated; the quoted value actually refers to 
the total excess of heavy-element material, above the so- 
lar abundance, in the entire planet. Erosion of the core 
and possible mixing of some core material into the enve- 
lope is not considered. This process has been shown to 
be unlikely for the case of Jupiter (Lissauer & Stevenson 
2007), but such estimates have not been extended to the 
case of planets in the 10 Mj up range. 

The structure of the hydrogen-helium envelope is cal- 
culated according to the differential equations of stel- 
lar structure (Kippenhahn & Weigert 1990), which as- 
sume hydrostatic equilibrium, a spherically symmetric 
mass distribution, radiative or convective energy trans- 
port, and energy conservation. The energy sources are 
provided by planetesimal accretion, contraction of the 
gaseous envelope, and cooling. The additional energy 
source provided by deuterium burning is included in the 
later phases of accretion and during the constant-mass 
final cooling phase, once the mass has exceeded 10 Mj up 
and internal temperatures exceed w 10 5 K. The full set 
of equations, supplemented by calculation of the mass 
accretion rates onto the core and the envelope, and of 
the planetesimal trajectories, is solved by the Henyey 
method (Henyey et al. 1964). 

At the inner boundary of the envelope the radius is 
set to i?corc, which is determined from its mean den- 
sity. During the earlier phases of the evolution, when 
the envelope mass is less than about 0.1 Mj up , the core 
is assumed to be a mixture of rock and ice with a mean 
density of 3.0gcm~ 3 . During the later phases, when the 
pressure at the base of the envelope increases to values 
above ~ 10 11 dynes cm -2 , an ANEOS equation of state 
with 50% rock and 50% ice (Thompson 1990) for the 
core is used to determine its mean density, which can 
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increase to 60 gem -3 or higher. In the hydrogen-helium 
envelope, the equation of state is taken to be given by the 
tables of Saumon, Chabrier, & van Horn (1995), which 
take into account the partial degeneracy of the electrons 
as well as non-ideal effects. The chemical composition is 
taken to be near-solar, with X = 0.70, Y = 0.283, and 
Z = 0.017, where X, Y, Z are, respectively, the mass 
fractions of hydrogen, helium, and heavy elements. The 
tables of course do not include a Z component, so the Y 
component was adjusted upwards to partially compen- 
sate. 

The Rosscland mean opacity during the formation 
phase combines the low-temperature atomic/molecular 
calculation of Alexander & Ferguson (1994) with the in- 
terstellar grain opacities of Pollack et al. (1985). The 
opacity values of the grain component are reduced by 
a factor 50 to approximately represent the reduction 
caused by grain growth and settling in the protoplanet 
(Podolak 2003; Movshovitz & Podolak 2008). However, 
in two of the runs the grain growth and settling are cal- 
culated in detail in the temperature range 100-1800 K 
as described in Movshovitz et al. (2010). The grain size 
distributions and the opacities are recalculated in every 
layer at every time step in that temperature range. These 
opacities are important in regulating the rate at which 
the envelope can contract, and therefore the rate at which 
it accretes gas. However, once the envelope is well into 
the rapid gas accretion phase, at about 0.25 Mj up , the 
gas accretion rate is limited by the physical properties of 
the protoplanetary disk near the planet, and the precise 
values of envelope opacity assume a less-important role. 
Once the planet reaches its final mass, say 12 Mj up , the 
grains are assumed to settle rapidly and to evaporate in 
the interior. For the final contraction/cooling phase at 
constant mass, the molecular opacities of Freedman et al. 
(2008) are used, with solar composition, up to a temper- 
ature of 3500 K. At and above that temperature, with 
any reasonable opacity, the interior is convective. 

At the outer surface of the envelope, the mass addi- 
tion rate of gas, during the earlier phases of accretion, 
is determined by the requirement that the planet radius 
R p match the effective accretion radius, which is given 
by (Lissauer et al. 2009) 



-Rcff 



GM P 
KR H 



(2) 



where c s is the sound speed in the disk, Rh is the 
Hill sphere radius, and M p is the total mass of the 
planet. The constant K 0.25 is determined by three- 
dimensional numerical simulations which calculate the 
accretion rate of gas from the protoplanetary disk onto 
the planet (Lissauer et al. 2009). As a result, in the limit 
where Rh is small compared with the Bondi accretion 
radius GM p /c 2 s , R cS = 0.25R H - 

Additional boundary conditions at the surface depend 
on the evolutionary phase. During the early phases when 
M p < 0.25Mj up , the density and temperature are set to 
constant values appropriate for the protoplanetary disk, 
Pnob and T ne |-,, respectively. The density p nc ^ is deter- 
mined from the assumed value of a using a standard 
gas-to-solid ratio of 70 and H p /a p = 0.05, where H p is 
the (gaussian) disk scale height and a p is the distance 
of the planet from the star. However at some point 



during the rapid gas accretion phase, the mass addition 
rate required by condition (2) exceeds the rate at which 
matter can be supplied by the disk. The disk-limited 
rates, based on three-dimensional hydrodynamic simu- 
lations, are described in the next section. During that 
phase, the boundary conditions at the actual surface of 
the planet, whose radius falls well below i? e ff , are deter- 
mined through the properties of the accretion shock at 
this surface, as described in detail by Bodenheimer et al. 
(2000). The basic assumption is that practically all of 
the gravitational energy released by the infalling gas is 
radiated away at the shock; this energy escapes through 
the infalling envelope ahead of the shock. This assump- 
tion defines the 'cold start' for planetary evolution. 

During the final phase of cooling at constant mass, the 
planet becomes isolated from the disk and the surface 
boundary conditions change again, to those of a black- 
body in hydrostatic equilibrium 



L 
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and 



kP = — q , 
3 y ' 



(3) 



where erg is the Stcfan-Boltzmann constant, T e g- is the 
surface temperature, L is the total luminosity, and re, 
P, and g are, respectively, the photospheric values of 
Rosscland mean opacity, pressure, and acceleration of 
gravity. Insolation from the star is not included. 

Significant deuterium burning in the mass range con- 
sidered begins near the end of the phase of rapid gas 
accretion. The burning occurs via the reaction 



2 D+ *H 



3 Hc + 7 



(4) 



with an energy release Qd P = 5.494 MeV per reaction. 
The initial deuterium abundance by mass fraction is set 
to 4 x 10~ 5 , consistent with the value derived from the 
local interstellar medium (Prodanovic et al. 2010). The 
reaction rate (reactions per second per gram) is taken 
from the Nuclear Astrophysics Compilation of Reaction 
Rates (Angulo et al. 1999): 



Rd P (p, Tq) = 



5.365 x lti 2% P X 1H X 2 H 



exp(-37.21/T 6 1/3 ) 



T. 



2/3 



{1 + T 6 [0.0143 + r 6 (3.95 x 10"'T e - 9.05 x 10 _t> )]},(5) 

where Tq is the temperature in 10 6 K, p is the density 
in cgs, Xih is the mass fraction of X H, and X 2 h is the 
mass fraction of 2 H (deuterium) . This rate is then mul- 
tiplied by the screening factor, which takes into account 
ion-ion and ion-electron screening in partially degener- 
ate dense plasmas (Potekhin & Chabrier 2012). The 
energy generation e, per gram per second, is then ob- 
tained, zone by zone, from the rate multiplied by Qdp 
in the appropriate units. To get the change in the deu- 
terium abundance during one time step, it is assumed 
that the planet interior is fully convective and therefore 
fully mixed. This assumption is valid for the planets 
considered during the phase of contraction and cooling, 
even if no deuterium is burned. The convective veloci- 
ties of order 10-100 cms -1 , calculated according to the 
mixing-length approximation, give a mixing time scale 
far shorter than the D-burning time scale. The reaction 
rate multiplied by zone mass is integrated over the entire 
envelope and used to calculate the abundance change. 
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Given the central stellar mass M*, the solid surface 
density a, and the distance of the planet from the star 
a p , the isolation mass for the solid material is 

M iso = ^Cf/Hl- 1/2 a^4 (6) 

where C is the number of Hill-sphere radii on each side 
of the planetary core from which it is able to capture 
planetesimals; C = 4 in our simulations. Once the core 
mass approaches M- mo , the dM core /dt slows down dras- 
tically, and beyond that point, gas accretion continues 
and surpasses the core accretion rate. The core mass in- 
creases to a value of about \^2Mi so at crossover, when 
-Mcore = Mcnv (Pollack et al. 1996). This phase of rela- 
tively slow accretion rates onto both core and envelope 
is known as 'Phase 2'. 

3. DISK-LIMITED GAS ACCRETION RATES 

The epoch of rapid gas accretion in the core-nucleated 
accretion model generally begins soon after the envelope 
mass, M env , exceeds the core mass, M core , as can also be 
shown by means of simple thcrmodynamical arguments 
(DAngclo et al. 2011). In a proto-solar nebula at ~ 
5 AU, this condition typically occurs when the planet 
mass M p = M core + M env is between ~ 10 to a few tens 
of Earth masses. After this point, the planet's envelope 
tends to contract very rapidly, limited only by the rate 
of energy escape at the surface, and a high rate of gas 
accretion is required to maintain the condition R p = i? e ff- 
At or about 0.25 Mj up this condition can no longer be 
met, the rate is set by the ability of the protoplanetary 
disk to deliver gas to the planet, and R p contracts well 
within i? ff. 

There are various regimes of disk-limited gas accretion 
(see DAngelo & Lubow 2008). For the purpose of this 
study, we are mainly interested in the high-mass limit 
Rh > Hp, where Rh = a p \J M v j (3A/+) is the Hill radius 
of the planet and H p is the disk thickness at the planet's 
orbital radius, a p . In this regime, disk- limited accretion 
rates can be affected by disk-planet gravitational interac- 
tions if tidal torques overcome viscous torques. Assume 
that the turbulent (kinematic) viscosity of the disk at the 
orbital distance of the planet is given by v t = atH p fl p , 
where £l p is the local Kcplcrian rotation frequency of 
the disk and at is the viscosity parameter. Then tidal 
torques exerted by the planet on the disk exceed viscous 
torques exerted by adjacent disk rings on each other if 

where A = max (H p , Rh) and / is a factor of order unity 
(see, e.g., DAngelo et al. 2011, and references therein). 
When the left-hand side of Equation (7) is much greater 
than the right-hand side, a gap forms in the disk surface 
density along the planet's orbital radius. 

We estimated disk-limited accretion rates, M p , us- 
ing high resolution 3D hydrodynamical simulations of 
a planet embedded in a protoplanetary disk. We used 
an approach along the lines of DAngelo et al. (2003). 
We considered a disk with a constant aspect ratio of 
H p /a p = 0.05 and with the parameter a t ranging from 
4 x 10~ 4 to 2 x 10~ 2 . The unperturbed surface density 
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Figure 1. Averaged surface density of gas in a disk, as a function 
of distance from the star, that tidally interacts with a planet with 
M p = 10 -2 M+. The orbital radius of the planet is a p . The surface 
density is plotted in scaled units of M* / a p . The aspect ratio of 
the disk is H p /a p = 0.05 and the turbulent viscosity parameter is 
a t = 10- 2 . 

of the disk is taken to be a power-law of the distance 
from the star with exponent —1/2. The planet is kept 
on a fixed circular orbit and the continuity and Navier- 
Stokes equations (written in terms of linear and angular 
momenta) for the gas are solved in a reference frame co- 
rotating with the planet. 

The disk is assumed to be vertically isothermal and, 
radially, the temperature drops as the inverse of the dis- 
tance from the star. At the radial distance a p , the tem- 
perature is T p = (/XdWH/^B^p^p; which is equal to 
53.8/^ K at 5AU from a solar-mass star (/j,d indicates 
the mean molecular weight of the disk's gas) . The rela- 
tively simple equation of state adopted here (the pressure 
p cx Tp, where p is the mass density and the tempera- 
ture T is a given function of the orbital distance) al- 
lows us to write the fluid equations in a non-dimensional 
form so that the gas accretion rates can be expressed in 
terms of a p Y. p fl p , where S p is the unperturbed gas sur- 
face density of the disk at the planet location (i.e., that 
the disk would have in the absence of the planet). Fur- 
thermore, the planet's mass enters the calculations only 
via its ratio to the mass of the star. We considered val- 
ues of the ratio M p /M* up to 0.02. The planet's gas 
accretion rate starts to decline for planet masses greater 
than the value for which the inequality in Equation (7) 
is satisfied. This critical mass is larger within more vis- 
cous disks. The equation also suggests that there is a 
dependence on the disk thickness, which, however, was 
not explored here. We notice that reasonable values of 
Hp/ a p for evolved disks, between 1 and 10 AU, range 
from ps 0.03 to ps 0.05 (e.g., DAngelo & Marzari 2012), 
affecting the right-hand side of Equation (7) by a fac- 
tor of less than 3 (for Rh > H p ), whereas uncertainties 
on at are much larger, spanning 2 orders of magnitude 
or more. An unperturbed surface density with a power 
index different from that adopted here (—1/2) may also 
affect the accretion rates. We expect these effects to be 
small, especially when tidal torques exerted by the planet 
drastically modify the surface density, which is typically 
the case in the models discussed here. 
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Figure 2. Disk-limited gas accretion rates as a function of the 
planet-to-star mass ratio, M p /M*, for H p /a p = 0.05 and for two 
values of the disk turbulent viscosity parameter: at = 4 X 10 
and 10 -2 . Symbols are data obtained from 3D hydrodynamics 
simulations: filled pentagons/open circles refer to the lower/higher 
viscosity case. The solid and dashed lines represent results from 
the fitting procedure outlined in the text. In the units of M p , T, p 
represents the unperturbed disk gas surface density and Q p the 
Keplerian rotation rate at the planet's orbital radius, a p . 

In the calculations, the disk domain extends in radius 
as close to the star as 0.1 a p (and 0.05 a p , in some calcula- 
tions) and as far as 9.4 a p . More vigorous perturbations 
exercised by larger mass planets cause the inner/outer 
disk radius to decrease/increase with increasing planet- 
to-star mass ratio. Figure 1 shows the surface density, 
averaged in azimuth, for a case in which M p = 10~ 2 M* 
and at = 10~ 2 . Notice that the low densities in the disk 
inside the orbit of the planet are a consequence of tidal 
torques and planet accretion (see Lubow & D'Angclo 
2006) , with possibly some impact from the finite radius of 
the grid inner boundary (0.05et p in the calculation shown 
in the figure). The analysis of Lubow & D'Angelo (2006), 
where applicable, suggests that the effects of the finite 
inner grid radius are small. 

High resolution in and around the planet's Hill sphere 
is achieved by means of multiple nested grids (DAngelo 
et al. 2002, 2003) centered on the planet's position. This 
methodology allows us to solve the fluid equations (lo- 
cally around the planet), and hence to resolve the accre- 
tion flow, on length scales of order 0.01 Rh, or w 7 Rj up 
at 5 AU. 

The gas that orbits the planet deep within its gravita- 
tional potential is eventually accreted. We assume that 
gas can be accreted within a spherical region of radius 
0.1 Rh (or 0.05 Rh in some models), centered on the 
planet. The amount of accreted gas is proportional to 
the amount of gas available in the region (see DAngelo 
& Lubow 2008, and references therein). In these calcu- 
lations, accreted gas is removed from the computational 
domain but not added to the mass of the planet in order 
to achieve a stationary accretion flow (see Lissauer et al. 
2009). 

We determined an interpolation procedure for the disk- 
limited gas accretion rates obtained from calculations, 
by performing piece-wise parabolic fits (in a logarith- 
mic plane) to each (M p , M p /M+) data set, relative to 
a given value of the turbulent viscosity. Two such fitting 



curves are shown in Figure 2 (explicit expressions are 
provided in Appendix A). Linear interpolations among 
these curves provide the accretion rate at the desired 
viscosity parameter, at- In doing so, we derived a func- 
tion M p = Mp(M p , M+, a p , T, p , a t ), which we employ in 
our planet formation calculations. We recall that E p here 
represents the disk gas surface density at the planet's or- 
bital radius, in the absence of the planet. An analytic 
formula is available for the accretion rate at the low-mass 
end (D'Angelo & Lubow 2008). However, in the forma- 
tion calculations the use of these curves is not required 
until M p exceeds ps 0.25 Mj up . 

During the disk-limited gas accretion phase, the solid 
accretion rate is arbitrarily limited to a fraction of the 
value at crossover; the precise value has practically no 
effect on the results. We don't expect the core-accretion 
prescription to be valid at this stage, because most solids 
in the disk will not be in the form of planctesimals, and 
we do not have the capability to model giant impacts. As 
the planet reaches within 2% of the desired final mass (e. 
g. 12 Mj up ) the gas accretion rate, already quite low, is 
smoothly reduced to zero. 

4. CALCULATIONS AND RESULTS 

A recent paper on deuterium burning in objects formed 
through the core-accretion scenario (Molliere & Mor- 
dasini 2012) considered the basic case of a body forming 
at 5.2 AU in a disk around a 1 M star with a solid sur- 
face density of a = 10 gem -2 and T ne b = 150 K. Their 
study compared results obtained by varying the follow- 
ing parameters: initial entropy of the object after forma- 
tion (hot start vs. cold start), helium abundance, metal 
abundance, initial deuterium mass fraction, a, which de- 
termines the final planet core mass, and maximum gas 
accretion rate. Their calculations differ from ours in the 
phase of rapid gas accretion, when disk-limited rates ap- 
ply. They take that rate to be an arbitrary parame- 
ter, while we use the three-dimensional simulations men- 
tioned above (Section 3) to determine it. Here we con- 
centrate on cold-start models and consider a somewhat 
different set of parameters: stellar mass, formation po- 
sition of the planet in the disk, solid surface density a, 
method of computation of the opacity in the planetary 
envelope during the formation phase, and protoplanetary 
disk viscosity parameter at- The planet's core mass is de- 
termined through the calculation itself, and it depends 
on the first three of these quantities. Note that the final 
core masses found in our calculations fall in the range 
4.8-31 M e , while those of Molliere & Mordasini (2012) 
are higher (30-100 M$). The formation and evolution 
are assumed to take place at a fixed orbital radius. 

The parameters for the runs are given in Table 1. The 
columns in the table give, respectively, the run identi- 
fier, the mass of the central star in M©, the distance of 
the planet from the star, the solid surface density a, the 
density p ne h at the surface of the planet during the ear- 
lier phases when this surface connects with the disk, the 
temperature T no b at the surface during the same phases, 
the method of opacity calculation during the formation 
phase — that is, whether it includes the calculation of 
grain settling and coagulation (gs) or not (ngs) — , the 
value of the viscosity parameter a t in the disk during 
the phases of disk-limited gas accretion, and the isola- 
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Table 1 

Input Parameters 



Run 


M/M 


Distance (AU) 


a (gem 2 ) 


Pneb (gem" 3 ) 


T neb (K) 


opacity 


at 




M iso (M e ) 


1A 


1 


5.2 


10 


9 x 10~ n 


115 


gs 


1.0 x 10~ 


-•> 


11.6 


IB 


1 


5.2 


10 


9 x 10~ n 


115 


ngs 


1.0 x 10~ 


-■> 


11.6 


1C 


1 


5.2 


4 


3.7 x 10~" 


115 


gs 


1.0 x 10~ 


-2 


2.9 


2A 


2 


9.5 


4 


1.8 x 10~ u 


125 


ngs 


1.0 x 10~ 


-•> 


12.6 


2B 


2 


9.5 


4 


1.8 x 10~" 


125 


ngs 


4.0 x 10~ 


-3 


12.6 


2C 


2 


9.5 


6 


2.8 x 10~ u 


125 


ngs 


1.0 x 10~ 


-•> 


23.2 



tion mass (Equation 6). 

Some results for the six runs are presented in Table 2. 
Each run is given two lines, the first for a final planet 
mass that burns less than half of its deuterium, the sec- 
ond for a nearby mass that burns more than half. The 
columns give, respectively, the run identification, the fi- 
nal planet mass in Mj up , the total time to reach the 
final mass (the formation time), the final core mass, the 
central temperature (T c .f, at the core/envelope interface) 
just after formation, the maximum central temperature 
during D-burning, the central density p c _f just after for- 
mation, the planet's radiated luminosity L c ^ just after 
formation, and the mass fraction of deuterium that re- 
mains after 4 Gyr of evolution, in units of the initial D 
mass fraction of 4 x 10 -5 . 

4.1. Results for 1 M Q 

Run 1A, with standard parameters of 1 M Q , 5.2 AU, 
and a = 10gcm~ 2 was originally calculated by 
Movshovitz et al. (2010) through most of the forma- 
tion phase, including the detailed calculation of grain 
opacity (their run crlO). Their run, whose characteristics 
are listed in that paper, ended at the beginning of disk- 
limited gas accretion, with a core mass of 16.8 M ffi and 
an envelope mass of 56.8 M ffi at a total elapsed time of 
1 Myr. In this work, the run was continued through the 
disk-limited phase with a t = 10 -2 (Section 3) up to the 
mass range required for deuterium burning. The max- 
imum gas accretion rate was 2.5 x 10 _1 M® yr" 1 at a 
total planet mass of 96 M ffi , declining to 10 -2 M e yr _1 
at 10 Mj up . For several different masses in that range, 
the accretion was terminated, the opacity was reset in 
the surface layers to the values given by Freedman et al. 
(2008), and the evolution was followed at constant mass 
up to Gyr times. The runs were terminated when deu- 
terium burning ceased, and the mass M50, where 50% 
of the original deuterium had been burned, was deter- 
mined. In the Run 1A, the total formation time at M50, 
up to termination of accretion, was 1.2 Myr, well within 
the lifetime of protoplanetary disks. 

The planetary luminosity as a function of time for three 
different final masses in Run 1A is shown in Figure 3, 
where it is compared with that typically obtained in a 
'hot-start' model. In the case of 16 Mj up , just after for- 
mation the central temperature T c j = 2.8 x 10 5 K, too 
low for substantial burning on a short time scale, even 
though p c ,i = 80gcm~ 3 . Under these conditions the 
screening correction factor to the nuclear reaction rate 
is high, about 88. Consequently, deuterium burning can 
take place at relatively low temperatures compared to 
those (~ 10 6 K) where deuterium burns in solar-mass 
stars. The central temperature T c as a function of time 
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Figure 3. Luminosity (in solar units) as a function of time for 
Run 1A during the post-formation deuterium-burning phase for 
three different planet masses. Solid curve: 16 Mj up , dashed curve: 
12.5 Mj U p, short-dash dot curve: 12 Mj up . The long-dash dot 
curve shows the results for a hot-start model of 10 Mj up (Baraffc 
et al. 2003). 

(at the core/envelope interface) gradually increases as a 
result of slow deuterium burning and is accompanied by 
a slight increase in radius. When T c reaches 3.2 x 10 5 K, 
a rapid increase in burning occurs, leading to a peak in 
luminosity at about 10 8 years. At the peak about 60% 
of the deuterium has burned, and T c is near its maxi- 
mum of 5.1 x 10 5 K. At the same time the radius has 
increased from 7.4 x 10 9 cm to 1.25 x 10 10 cm; then it 
contracts again after the luminosity peak. At the end of 
the evolution essentially all the deuterium has burned. A 
similar process, involving a rapid increase in deuterium 
burning in the context of a slowly accreting brown dwarf, 
was studied by Salpeter (1992); he denotes the event a 
'deuterium flash'. The radii for the three masses, as well 
as for the hot-start case, are shown in Figure 4. The 
general result that cold-start models result in a radius 
increase during deuterium burning agrees with the pre- 
vious results of Molliere & Mordasini (2012). 

In the case of 12.5 Mj up , right after formation the cen- 
tral temperature is lower, only 2.6 x 10 5 K, with a central 
density of 52 g cm -3 and a screening factor of 70. It takes 
almost 10 9 years for rapid deuterium burning to start, at 
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Figure 4. Radius (in 10 10 cm) as a function of time for Run 1A 
during the post-formation deuterium-burning phase for three differ- 
ent planet masses. Solid curve: 16 Mj up , dashed curve: 12.5 M,j up , 
short-dash dot curve: 12 Mj up . The long-dash dot curve shows the 
results for a hot-start model of 10 Mj up (Baraffc et al. 2003). 



T c = 2.8 x 10 5 K, with the burning occurring on a much 
longer time scale than in the case of 16 Mj up . Eventually 
about 98% of the deuterium is burned, and the luminos- 
ity peak, which is somewhat lower, is shifted to later 
times. At the peak, about half of the D has burned, and 
this point is also close to the maxima in radius and T c . In 
the case of 12 Mj up only 6% of the deuterium is burned, 
and no peak in luminosity appears. At M50 itself, the 
peak involves only a factor 2 increase in luminosity. In 
the peaks, the total energy J Ldt is found to agree closely 
with the total energy available from D-burning, given by 
the quantity Qdp/nid X Mp x Xdfd, where Qdp, the en- 
ergy production per reaction is expressed in ergs, is 
the mass of a deuterium atom, M p is the planet mass, 
Xd is the initial mass fraction of deuterium, and fd is the 
fraction of the initial D that burned. The total energy is 
about 2.5 x 10 45 ergs for the 12.5 Mj up case. 

Results for two runs whose final masses closely bracket 
M50 are shown in Table 2. The main result of this case 
is that M50 = 12.37 Mj up with a heavy-element core 
mass of 16.8 M©. By way of comparison, a cold start 
model with a core, calculated by Molliere & Mordasini 
(2012) with about the same basic parameters (1 M©, 
a = lOgcm , a p = 5.2 AU), with a similar helium mass 
fraction of 28%, but with some differences in assumptions 
and computational procedure, gives M50 = 12.6 Mj up . 

The maximum T c at the core/envelope interface for 
M50 in this case is close to 3.2 x 10 5 K, a very sensitive 
function of mass. Whether significant D-burning occurs 
depends sensitively on this temperature. If it reaches, say 
2.5 x 10 5 K, practically no D is burned for the correspond- 
ing mass of 12.0 Mj up . If it reaches 4.0 x 10 5 K, practi- 
cally all (98%) of the D is burned for the corresponding 
mass of 12.5 Mj up . Once the threshold is reached, en- 



ergy deposition from burning increases the temperature, 
which increases the reaction rate, as it is proportional 
to T 12 . The resulting expansion leads to a near thermal 
equilibrium, with the energy produced from D-burning 
matched closely by the total radiated luminosity. 

Run IB differs from 1A only with respect to the cal- 
culation of the opacity resulting from grains in the pro- 
toplanetary envelope during the formation phase. As 
mentioned above, in Run 1A this opacity is obtained 
through detailed consideration of grain settling and co- 
agulation (Movshovitz et al. 2010). In IB a table of 
interstellar grain opacities is used, reduced by a factor 
of about 50. The characteristics of this run, up to a 
mass of about 1 Mj up , are very similar to those listed for 
Run lsG in Lissauer et al. (2009). The crossover mass 
is 16.16 M©, the crossover time is 2.31 Myr, and the on- 
set of disk-limited rapid gas accretion occurs at a core 
mass of 16.8 M© and a time of 2.41 Myr. Note that the 
evolution time up to this point is 2.4 times longer than 
in Run 1A. Note also that the core mass is the same as 
in Run 1A; the substantial difference in opacity, which 
can be up to two orders of magnitude in certain (/?, T) 
regions, has practically no effect on the core mass. 

Here, the disk-limited accretion rates are used to con- 
tinue the evolution up to the D-burning mass range. The 
luminosity as a function of time up to the end of accre- 
tion is shown in Figure 5. The results for D-burning after 
that time show that M 50 = 12.20 Mj up , not significantly 
different from the results of Run 1A. As Table 2 shows, 
T Cj f in Run IB, at the same final mass, is slightly higher 
than that in Run 1A, just after formation. Correspond- 
ingly, p C: f is slightly lower. These small differences indi- 
cate a slightly higher entropy for IB after formation, as 
indicated by the slightly higher luminosity at this point. 
The increased envelope opacity in IB as compared with 
1A results in slower heat loss and tends to keep inter- 
nal temperatures higher. However this effect is almost 
compensated by the fact that the formation time is more 
than twice as long in IB. Even the slight increase in T c .f 
in IB as compared with 1A allows M50 to be pushed to 
a slightly lower mass. 

Run 1C differs from Run 1A in that a is reduced 
to 4 gem -2 , a value only slightly greater than that in 
a minimum-mass solar nebula (Weidenschilling 1977). 
Grain settling and coagulation are included in the opac- 
ity calculation. The earlier portions of this run, up to 
the onset of disk-limited gas accretion, are described in 
Movshovitz et al. (2010), their run <r4. The time to reach 
this point, 3.5 Myr, is considerably longer than in Run 
1A, first, because the core accretion rate is considerably 
lower, and second, because the lower isolation mass re- 
sults in reduced luminosity and reduced gas accretion 
rate during Phase 2 (Pollack et al. 1996). The crossover 
mass is 4.09 M©, and the core mass at the time of onset 
of disk-limited accretion is 4.74 M©. 

The calculations were continued up to the point where 
gas accretion terminated, at which point the core mass 
was 4.8 M©. The total time to reach 12 Mj up was about 
4.1 Myr, and to 14 M Jup , about 4.5 Myr. The peak 
disk-limited accretion rate was 1.0 x 10 -1 M© yr -1 at 
0.3 Mj up , a factor of 2.5 lower than in Run lA because 
of the reduction in S p by the same factor. By the time 
the total mass was 5 Mj up the rate was down to 1.8 x 
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10 2 M® yr 1 , and at 10 Mj up it had further declined to 
4.3 x 10 -3 M ffi yr -1 . Much of the time during the disk- 
limited accretion phase was spent in accreting the last 1— 
2 Mj up to reach the D-burning point. The luminosity as 
a function of time for this run, up to the end of accretion, 
is shown in Figure 5. 

The luminosity versus time plots for Run 1C during 
the D-burning phase look similar to those for 1A, except 
in this case M50 noticeably increases to 13.55 Mj up . The 
reduced core mass in 1C (4.8 M®) as compared to that 
in 1A (16.8 M®) is clearly associated with the difference, 
in agreement with the results of Molliere & Mordasini 
(2012). In our calculations, the core equation of state 
gives a core radius of 3.8 x 10 8 cm for the Run 1C core 
of mass 4.8 M® when the total mass is 12 Mj up . For the 
core of 16.8 M® in Run 1A, at the same total mass, the 
radius is 6.0 x 10 8 cm. Thus, at the core boundary, the 
gravitational potential is more negative, and the gravity 
is about 40% greater in 1A than in 1C. The calculated 
values of T Cif are w 2.6 x 10 5 K and w 2.1 x 10 5 K in 
Runs 1A and 1C, respectively. 

It follows from the equation of hydrostatic equilibrium 
(Molliere & Mordasini 2012) that in a convective enve- 
lope the adiabatic temperature gradient at the interface 
should be proportional to the core gravity, so a higher 
gravity most probably gives a higher temperature. How- 
ever, this statement is inconclusive. We calculated static 
models for a planet of 12 Mj up , all with the envelope 
entropy of Run 1C, with core masses ranging from to 
15 M®. We found practically no difference in T c as a 
function of M core , with T c decreasing by less than 1% 
when the core mass increases from to 15 M®. 

The real source of the difference in T c f between Runs 
1A and 1C is the entropy in the envelope. The lower T Cj f 
and higher p c j for 1C as compared with 1A indicate a 
lower entropy, which is consistent with the fact that the 
luminosity just after formation is lower by more than a 
factor 2 in 1C (Table 2). The values of entropy just after 
formation for a planet of 12 Mj up in Runs 1A and 1C are, 
respectively, 8.02 and 7.52 ks per baryon. The entropy 
is determined through the physical processes that occur 
during the entire formation phase; for example, the for- 
mation time for Run 1C is almost 4 times longer than 
that for 1A, and the same opacities were used, which 
suggests a lower entropy. Thus there exists a qualitative 
understanding of the relation between core mass and T Cj f , 
but a quantitative theory, apart from the numerical sim- 
ulations, is quite difficult. 

In Run 1C, T Cj f = 2.1 x 10 5 K is the maximum reached 
for a final mass of 12 Mj up , and it is insufficient for D- 
burning. In the case of Run 1A, the corresponding T c j is 
much closer to the threshold required for burning. Thus 
the planet with the higher M core is able to produce signif- 
icant D-burning at a lower total mass. As Table 2 shows, 
in the mass range for Run 1C where D-burning begins, 
just above 13.5 Mj up , T Cj f is somewhat less (2.3 x 10 5 K) 
than in the corresponding mass range for Run 1A. How- 
ever, to compensate, p c> f is higher, about 65 gem , and 
the screening factor at the center has increased to 160. 
Again, the lower entropy at formation for 1C, as com- 
pared with 1A, a result of various processes associated 
with the accretion of core and envelope, leads to a higher 
M 50 . 
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Figure 5. Luminosity (in solar units) as a function of time for 
Run IB (long- dashed curve), Run 1C (short- dash long- dash curve), 
Run 2A (solid curve), and Run 2C (short-dashed curve) during the 
formation phase. 

4.2. Results for 2 M 

The formation phases of Runs 2 A and 2 C, for a central 
star of 2 M®, are illustrated in Figure 5, which gives the 
luminosity as a function of time, and Figure 6, which 
gives the core mass, envelope mass, and total mass as 
a function of time. Run 2A differs from 1A in that the 
planet is placed 9.5 AU away from a star of 2 M®, in a 
disk with a — 4 g cm' 2 . In a minimum mass solar nebula, 
scaled to the mass of this star, the corresponding value 
would be 2gcm~ 2 . 

The isolation mass, however is quite similar to that in 

IA, 12.6 rather than 11.6 M®. The opacity during the 
formation phase of 2A is taken from a table of interstel- 
lar grain opacities, reduced by a factor of 50, as in Run 

IB. However, the comparison between 1A and IB showed 
that these opacities have little effect on A/50. The forma- 
tion time is longer in 2A than in 1A because of the longer 
dynamical time at the larger distance, the reduced solid 
surface density, and the somewhat higher envelope opac- 
ity. However, these effects are partially compensated for 
by the smaller planetesimal size (50 km in 2A; 100 km 
in 1A), which increases the capture cross section 7ri? 2 apt , 
and by the increased gravitational focussing factor F g at 
the larger distance. 

The first luminosity peak for Run 2A (Figure 5) oc- 
curs at t = 3.54 x 10 5 yr, with log L/L® = -5.14, 
with M coro = 9.3 M®, with M onv = 0.024 M®, and with 
Mcore = 6.67 x 10~ 5 M® yr -1 . This peak corresponds to 
the maximum in the accretion rate of solids onto the core. 
The crossover mass (Figure 6) of 17.6 M® is reached in 
2.0 x 10 6 years. The second, much higher luminosity peak 
at 2.07 x 10 6 yr corresponds to the phase of rapid gas ac- 
cretion up to a final mass of 15 Mj up . At that time the 
maximum gas accretion rate is 2.2 x 10 _1 M® yr -1 and 
M p = 0.47 Mj up . Formation is complete, up to 15 Mj up , 
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Table 2 

Selected Results 



Run 


Af flna l/Mj u p 


iform (yr) 


M c ore/Me 


%,f (K) 


T max (K) 


Pel (gem -3 ) 


log (L f /L ) 


-Dfinal/Amt 


1A 


12.26 


1.19 x 10 e 


16.8 


2.60 x 10 5 


2.60 x 10 5 


49.7 


-6.22 


0.895 


1A 


12.48 


1.20 x 10 6 


16.8 


2.62 x 10 5 


3.54 x 10 5 


51.4 


-6.23 


0.164 


IB 


12.14 


2.67 x 10 e 


16.8 


2.69 x 10 5 


2.69 x 10 5 


47.0 


-5.84 


0.860 


IB 


12.26 


2.68 x 10 6 


16.8 


2.72 x 10 5 


3.18 x 10 5 


47.7 


-5.91 


0.328 


1C 


13.5 


4.37 x 10 6 


4.83 


2.31 x 10 5 


2.31 x 10 5 


65.1 


-6.57 


0.938 


1C 


13.6 


4.39 x 10 6 


4.83 


2.33 x 10 5 


4.27 x 10 5 


65.9 


-6.57 


0.292 


2A 


11.9 


2.14 x 10 e 


18.7 


2.73 x 10 5 


2.81 x 10 5 


41.4 


-5.64 


0.767 


2A 


12.0 


2.14 x 10 6 


18.7 


2.76 x 10 5 


3.06 x 10 5 


41.8 


-5.63 


0.435 


2B 


12.0 


3.23 x 10 e 


18.8 


2.78 x 10 5 


2.87 x 10 5 


44.6 


-5.73 


0.582 


2B 


12.1 


3.26 x 10 e 


18.8 


2.80 x 10 5 


3.28 x 10 5 


45.1 


-5.72 


0.320 


2C 


11.6 


8.75 x 10 5 


31.0 


3.37 x 10 5 


3.46 x 10 5 


27.9 


-4.80 


0.670 


2C 


11.7 


8.75 x 10 5 


31.0 


3.39 x 10 5 


3.56 x 10 5 


28.3 


-4.80 


0.390 
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Figure 6. Planet mass (in Mg)) as a function of time for Run 
2A and Run 2C during the formation phase. The final masses are 
15 Mj up and 12 Mj up , respectively. For Run 2A, the solid curve 
gives the core mass, the dotted curve the envelope mass, and the 
short-dash dot curve the total mass. For Run 2C, the short-dashed 
curve gives the core mass, the long-dashed curve the envelope mass, 
and the long-dash dot curve the total mass. 

in 2.2 x 10 6 yr. The final M core = 18.7 M e is slightly 
higher than in Run 1A. 

The luminosity as a function of time during the later 
deuterium-burning phase is shown for five different fi- 
nal masses in Run 2A in Figure 7. In the cases of 16 
and 15 Mj up , practically all (> 99%) of the deuterium is 
burned; in the case of 13.5 Mj up , about 92% is burned; 
for 13.0 Mj up , 75% is burned, and for 12.0 Mj up , just over 
50% is burned. Thus the value for M 50 « 11.95 M Jup is 
very close to the values obtained for Runs 1A/1B despite 
substantial differences in assumptions and initial condi- 
tions. As discussed in the comparison between Runs 1A 
and 1C, the somewhat larger M COTe in 2 A as compared 
to 1A is the main reason for the slightly lower M50 in 
2A. After formation, 2A has a slightly higher T Cj f than 
1A and slightly lower p Cj f, leading to a slightly higher 
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Figure 7. Luminosity (in solar units) as a function of time for 
Run 2A during the post-formation deuterium-burning phase for 
five different planet masses. Long-dash dot curve: 16 Mj up , solid 
curve: 15 Mj up , short-dashed curve: 13.5 Mj up , long-dashed curve: 
13 Mj up , short-dash dot curve: 12 M,j up . 



entropy. 

Comparing the luminosity curves for a mass of 16 Mj up 
in Figures 3 and 7, they look very different but in fact 
they are consistent. In Run 2 A (Figure 7) the higher 
T Ci f (because of the somewhat higher M core ) allows D- 
burning to start earlier than in Run 1A, and the value of 
L at the starting point is a factor of 4 higher. In fact the 
full widths of the two curves are quite similar, the peak 
values agree to better than a factor 2, and the integrated 
luminosities over time of the two curves agree to within 
10%. 

Run 2B has exactly the same parameters as 2A except 
that at is reduced by a factor 2.5, which affects the gas 
accretion rates during the disk-limited phase. Thus the 
formation time in 2B turns out to be a factor of 1.5 longer 
at 3.2 x 10 6 years, but still within the range of observed 
disk lifetimes. Table 2 shows that for final mass 12 Mj up 
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Figure 8. Luminosity (in solar units) as a function of time for 
Run 2C during the post-formation deuterium-burning phase for 
four different planet masses. Solid curve: 15 Mj up , long-dashed 
curve: 13.7 Mj up , short-dash dot curve: 12 Mj up , short-dash 
curve: 11.7 Mj up . The long-dash dot curve shows a hot-start model 
for 10 Mj up (Baraffc ct al. 2003). The cross gives the position and 
error bars for the companion to Beta Pic (Bonnefoy et al. 2013). 

the fractions of deuterium burned are in agreement for 
runs 2 A and 2B, within the uncertainties of the calcula- 
tions. Thus at has practically no effect on M50 in this 
case. Run 2B has a slightly lower entropy than 2A at 
12 Mj up , 8.2 fee per baryon versus 8.25, and therefore 
a slightly higher M50. Thus it appears that the longer 
time during disk-limited accretion in Run 2B has only a 
weak effect on both M50 and the entropy, at the same 

-^core • 

Run 2C has the same parameters as 2 A except that 
the solid surface density a is increased by a factor 1.5 to 
6 gem -2 . The first luminosity peak (Figure 5) occurs at 
t = 3.07 x 10 5 years with log L/L e = -4.45 at M core = 
15.6 M®. The crossover mass (Figure 6) is reached at 
t = 7.88 x 10 5 yr with a value of 30.7 M e . The maximum 
luminosity in the second peak is above log L/L Q = —1, at 
t = 7.915 x 10 5 years and a total mass of 0.62 M Jup . The 
higher a with respect to Run 2A results in a markedly 
higher M corc = 31 M® and a markedly shorter formation 
time (8.75 x 10 5 yr at M50). Despite these relatively large 
differences, the value for M 50 in 2C is only 2.5% smaller 
than in 2A. At the end of the formation phase, central 
temperatures are higher and central densities are lower 
in 2C as compared with 2A. Also, the screening factor is 
only 14 in 2C compared with 41 for 2A. The entropy at 
formation, for a final mass of 12 Mj up , is higher in 2C, 
9.08 fcs/baryon as compared with 8.25, corresponding 
to a higher luminosity at that point. The slope in the 
(M C ore, M 50 ) diagram between M cme = 18.7 and 31 M® 
is -0.024, a result which differs somewhat from that of 
Molliere & Mordasini (2012). They obtain a slope (in 
the same units) of -0.01, although for a different core 
mass range, 30 to 100 M®. 



Plots of luminosity versus time are shown in Figure 8 
for several different masses in Run 2C. As in Figure 7 
the higher masses give higher peak luminosity at earlier 
times than the lower masses, and at M50 there is only a 
very small peak. The L(t) curve for 15 Mj up starts at 
a higher value and reaches a maximum sooner than for 
the same mass in Run 2A, because of the higher internal 
temperature, but the value of log L at the peak is about 
the same. At 13.7 and 15 Mj up practically all of the 
initial D is burned. At 12 Mj up , 72% is burned, while at 
11.7 Mj up , which is very close to M 50 , 61% is burned. 

Plots of T c versus time, during the deuterium-burning 
phases, are shown for masses 12 and 15 Mj up in Figure 9, 
where they are compared with the results from Run 2A. 
The plot shows the effect of varying the core mass at 
fixed total mass, and of varying the total mass at fixed 
core mass. For example, for Run 2A at 15 Mj up the 
maximum T c is 4.6 x 10 5 K, while for 12 Mj up it is only 
3.06 x 10 5 K and is reached at a much later time. The 
vertical portions of these curves show the effect of rapid 
gas accretion from about 1 Mj up to the final mass of 
either 12 or 15 Mj up . The two nearly horizontal curves 
are for a total mass of 12 Mj up and core masses of 18.7 
(lower; Run 2A) and 31 M e (upper; Run 2C). The higher 
core mass results in a higher temperature by a factor of 
about 1.3. In the case of the 31 M® core, about 75% 
of the deuterium is burned; in the 18.7 M® core, a little 
over 50%. Note that the D-burning occurs late in the 
evolution, where small peaks in the temperature are seen. 
The remaining two curves correspond to a total mass of 
15 Mj up , with the same two core masses just mentioned. 
The D-burning occurs earlier than in the case of the lower 
total mass, and the higher core mass again gives a higher 
maximum T c . In both cases for 15 Mj up practically all 
the D is burned, and the residual mass fraction is smaller 
(8.9 x 10 -11 ) for the higher core mass as compared to the 
lower (2.6 x 10~ 9 ). 

In Run 2C with 15 M Jup , T c goes up to about 5 x 10 5 K 
and there are actually two minor peaks. Figure 10 illus- 
trates in more detail how various quantities vary during 
this phase. In this case, with a high T c .t , nuclear burning 
starts very early. During most of the phase, the object 
is not in thermal equilibrium. The first maximum in T c 
occurs when about 25% of the D has burned, close to the 
time of the maximum in the nuclear burning luminosity 
L nuc . Here L n uc is well above the radiated luminosity 
L, and the extra power goes into expansion, resulting in 
slight cooling of the interior. When half the deuterium 
has burned (1.3 x 10 7 yr), there is a maximum in lumi- 
nosity and radius, corresponding to the slight minimum 
in T c . Then contraction along with a slow decrease in 
nuclear burning leads to slight heating, and the second 
maximum occurs when 98% of the D has burned. This 
maximum corresponds to the time when L nuc starts to 
drop rapidly and to fall well below L. Beyond that point, 
even though contraction is occurring, there is insufficient 
burning to maintain the high temperature, and the ob- 
ject enters its final cooling phase. In contrast, in the case 
of 12 Mj up , the main D-burning in Run 2C takes place at 
practically constant T c , radius, and L, with a slight max- 
imum in T c of 3.67 x 10 5 K at about 10 8 yr. In this case 
the configuration is close to thermal equilibrium through 
most of the D-burning phase. 
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Figure 9. Central temperature T c (at the core/envelope inter- 
face) as a function of time for four cases. Solid curve: Run 2A 
(15 Mj up ), long-dashed curve: Run 2A (12 Mj up ), dash-dot curve: 
Run 2C (15 M Jup ), short-dashed curve: Run 2C (12 Mj up ). The 
phases of rapid gas accretion and final phases of constant mass 
with deuterium burning are shown. Core masses for Runs 2A and 
2C are 18.7 and 31 M®, respectively. 
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Figure 10. Detail of the deuterium-burning phase for Run 2C, 
15 Mj up . Solid curve: outer radius R p as a function of time, in 
units of 5 X 10 9 cm; long-dashed curve: central temperature T c , at 
the core/envelope interface, as a function of time, in units of 10 5 K; 
dash- dot curve: nuclear luminosity as a function of time in units of 
log (Lmic/L©) +7; short-dashed curve: radiated luminosity L as a 
function of time, in the same units as L nuc . 



4.3. Comparison with beta Pictoris b 

The cross in Figure 8 gives the approximate location 
of the directly imaged companion (Lagrange et al. 2009) 
to the well-known star Beta Pictoris. That star, ac- 
cording to http : / / exoplanet . eu, has a mass of about 
1.8 M and an age of 12 (+8,-4) Myr. The planet is 
located between 8 and 15 AU from the star (Lagrange 
et al. 2010); thus an approximate comparison can be 
made with these calculations. Beta Pic b's position in 
the (log L, t) diagram is plotted in Marois et al. (2010) 
where it is shown to fall on a theoretical track with mass 
10 Mj up as calculated from a 'hot start' by Baraffe et al. 
(2003). In http : //exoplanet . eu that mass is given as 
8 (+5,-2) Mj up . The surface temperature T e g has been 
estimated from observed near infrared colors (Bonncfoy 
et al. 2011; Quanz et al. 2010) at 1700 K, with consid- 
erable uncertainty (~ 300 K). Further infrared and as- 
trometric observations (Bonncfoy et al. 2013) are essen- 
tially in agreement, giving log (L/L Q ) = — 3.87± 0.08, 
Tes = 1700 ± 100 K, a p = 8 - 10 AU, and 'hot-start' 
masses in the range 7-13 Mj up . The bolometric luminos- 
ity found by Marleau & Cumming (2013) is in agreement 
with the above value, and they find 'hot-start' masses in 
the range 7-12 Mj up . 

In our 'cold-start' calculations the track for Run 2C, 
13.7 Mjup, passes close to the object in the (log L, t) 
diagram, and the calculations give T e g = 1627 K at an 
age of 12 Myr. Our mass 10 Mj up cannot possibly pro- 
vide a fit. The 'hot-start' models thus would show that 
the object is a planet, as defined by an object with mass 
not high enough to burn deuterium. However this par- 
ticular 'cold-start' model indicates that beta Pictoris b 
is presently burning deuterium, which, according to the 
same definition, would classify it as a brown dwarf. As 
mentioned in Section 1, this definition is not universally 
agreed upon; an alternative definition, based on the min- 
imum in the mass distribution of low-mass companions, 
observed within several AU of sunlikc stars, places the 
limit at « 25 Mj up . In this case Beta Pic b would still 
be a planet. Note that in the 'cold-start' calculations, 
the fit at 13.7 Mj up with an assumed a = 6 g cm" 2 is not 
unique; the companion could also be fit at a = 4 gem 
at a slightly higher mass, about 15.6 Mj up . Furthermore, 
these masses are uncertain and will probably change 
when the calculations are redone in the future with more 
detailed model atmospheres. Nevertheless, as such they 
are marginally consistent with the upper limits to the 
mass of Beta Pic b derived from radial velocity measure- 
ments (Lagrange et al. 2012). For a planet at 9 AU the 
limit is 12 M Jup ; at 10 AU it is 15.4 M Jup . 

We note also that the luminosity curve for 11.7 Mj up 
in Figure 8 agrees well with the observed luminosity of 
the directly imaged planet HR 8799 c at the stellar age 
(«6x 10 7 years). The observed value is given by Marlcy 
et al. (2012) as log L/L Q = -4.9 ± 0.1. The agree- 
ment of course requires a core mass of ps 30 M®. A 
hot-start model of about 10 Mj up without a core also 
agrees. However we do not make a detailed comparison 
with HR 8799 c, because the metallicity of the star is 
low ([Fe/H] = -0.47) and the planet orbits at 43 AU, 
making it highly debatable whether it could have formed 
by core-nucleated accretion. 
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Table 3 

Summary 



Run 


M/M Q 


Distance (AU) 


a (gem -2 ) 


Mcore (M ffi ) 


M 50 (M Jup ) 


1A 


1 


5.2 


10 


16.8 


12.37 


IB 


1 


5.2 


10 


16.8 


12.20 


1C 


1 


5.2 


4 


4.83 


13.55 


2A 


2 


9.5 


4 


18.7 


11.95 


2B 


2 


9.5 


4 


18.8 


12.05 


2C 


2 


9.5 


6 


31.0 


11.65 



5. SUMMARY AND CONCLUSIONS 

We investigate the boundary between brown dwarfs 
and giant planets, according to the definition that brown 
dwarfs can burn the deuterium that is present when they 
form, and giant planets cannot. The main parameters 
and the results for M50, the boundary mass at which 
half of the original deuterium is burned after 4 Gyr, are 
summarized in Table 3. The columns give, respectively, 
the run identification, the stellar mass in M©, a p , the 
initial disk solid surface density a at a p , the resulting 
M core , and M 50 . The main cases considered involve a 
planet/brown dwarf at 5.2 AU around a solar-mass star, 
and a planet/brown dwarf at 9.5 AU around a star of 
2 M©. The table shows that there is only a small varia- 
tion in the values of -M50, which, however, correlate with 
the core mass in the sense that the smaller the core mass, 
the higher the value of M50. 

The calculations, taken as a whole, indicate that the 
envelope entropy, which is a function of initial conditions 
and which is closely related to the core mass through 
the accretion processes during the formation phase, is an 
important factor in determining M50. However, certain 
physical processes during formation are shown to have 
only a small effect. Run IB has the same parameters as 
Run 1A except that the dust opacity during the forma- 
tion phase is higher by a factor that ranges from 2 to 100, 
depending on the depth in the envelope. This difference 
has a negligible effect on M50. Run 2B has the same pa- 
rameters as 2A except that the disk viscosity during the 
phase of disk-limited gas accretion is lower by a factor 
2.5. This difference also has a negligible effect on M50 
However the disk viscosity is important in another re- 
spect. If it is significantly lower than the range presented 
here (at ~ 10~ 2 ), then there will not be time to accrete 
a planet with mass necessary to burn deuterium during 
the lifetime of the disk. The gas accretion rate onto a 
planet of 4 Mj up around a star of 2 M© , in a disk with 
at — 4 x 10~ 4 , is reduced by a factor 400 compared with 
a disk with at = 4 x 10~ 3 (Lissauer et al. 2009), corre- 
sponding to less than a Jupiter mass in a million years for 
the initial conditions of Run 2B (formation time about 3 
Myr) . Of course the minimum viscosity required to build 
a planet up to about 12 Mj up will depend on parameters 
such as a and a p . This question is discussed in more 
detail in Appendix A. 

Core accretion models, in the cold-start case, are 
known to have low entropy compared with hot-start mod- 
els. In Marley et al. (2007) the entropy just after forma- 
tion for 10 Mj up was found to be 8.2 Ub per baryon for 
M cme = 16-8 M®. The corresponding luminosity at ages 
of 10 7 to 10 8 years was about 2 x 10 6 L©, certainly 




10 20 30 

Core Mass (M E ) 

Figure 11. The entropy in the interior of planets of total mass 
12 Mj up , immediately after formation, is plotted against their core 
masses, in Mq. The points plotted, from top to bottom, are from 
Runs 2C, 2A, IB, and 1C. 



fainter than observed values for directly imaged planets. 
In this mass range, for the given core mass, the entropy 
is very insensitive to the planet's total mass, as shown in 
that paper and confirmed by the present results. How- 
ever our calculations show that the entropy is quite sensi- 
tive to the core mass, as illustrated in Figure 11 (a similar 
effect has been found independently by Mordasini (2013) 
for M core > 20 M©). The points shown are all calculated 
with the same total mass and the same disk viscosity. 
All used the reduced interstellar grain opacity, except for 
the point at M core = 4.8 M©, for which the grain-settling 
opacities were used (if the interstellar opacities had been 
used, the formation time would have been considerably 
longer). However the comparison between Runs 1A and 
IB, which looked at the effect of changing the opacities, 
showed that the difference in entropy was less than 0.1 
ks per baryon at the same total mass. The effect on the 
entropy of changing the viscosity (Runs 2 A and 2B) was 
even smaller. Physical effects that do affect the entropy 
include the planetesimal accretion rate and the rate of 
contraction of the envelope, both of which affect the in- 
ternal heating of the envelope. Thus the luminosities of 
newly formed massive planets, depending on formation 
conditions, can vary by up to two orders of magnitude. 
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Table 4 

Data for Figure 11 



Run 


M fina i/M. Iup 


iform (yr) 


M 00 re/M® 


Tc,f (K) 


T max (K) 


Pec (gem -3 ) 


log (Lf/L ) 


-Dflnal/Anit 


2C 


12.0 


8.83 x 10 5 


31.0 


3.48 x 10 5 


3.67 x 10 5 


29.6 


-4.80 


0.278 


2A 


12.0 


2.14 x 10 e 


18.7 


2.76 x 10 5 


3.06 x 10 5 


41.8 


-5.63 


0.435 


IB 


12.0 


2.66 x 10 e 


16.8 


2.66 x 10 5 


2.66 x 10 5 


46.1 


-5.91 


0.930 


1C 


12.0 


4.10 x 10 6 


4.80 


2.08 x 10 5 


2.08 x 10 5 


54.0 


-6.56 


1.000 



Information on the runs whose entropies are plotted in 
the figure is given in Table 4. The table is in the same for- 
mat as Table 2 and gives the runs in order of decreasing 
entropy. Clearly, for this set of models, a lower entropy is 
associated with a longer formation time. The luminosity 
plots for these four cases in Figure 5 illustrate the same 
effect. 

The combination of M*, a p , and a determines the iso- 
lation mass, and thereby the ultimate core mass, which 
turns out to be a key factor in determining the entropy 
of the planet at formation. Higher entropy, in partic- 
ular the higher temperature, favors more rapid nuclear 
burning, so the higher entropy runs result in lower val- 
ues of M50. Nevertheless, the range of initial conditions 
explored here, which is considerable, produces only a 
small range in M50, about 11.6-13.6 Mj up , in agreement 
with previous independent calculations. We can further 
conclude, that for cold-start core-accretion models that 
do burn deuterium, the tracks in the luminosity versus 
time diagram can potentially provide agreement with the 



properties of directly-imaged low-mass stellar compan- 
ions. 

Primary funding for this project was provided by 
the NASA Origins of Solar Systems Program grant 
NNX11AK54G (P. B., G. D., J. L.). CD. acknowledges 
additional support from NASA grant NNX11AD20G. 
P. B. acknowledges additional support from NSF grant 
AST0908807. D. S. is supported in part by NASA grants 
NNH11AQ54I and NNH12AT89I. The authors are in- 
debted to Gilles Chabrier for the use of his nuclear 
screening factors. The 3D hydrodynamical simulations 
reported in this work were performed using resources 
provided by the NASA High-End Computing (HEC) 
Program through the NASA Advanced Supcrcomputing 
(NAS) Division at Ames Research Center. G.D. thanks 
Los Alamos National Laboratory for its hospitality. The 
authors thank the referee Dr Christoph Mordasini for a 
detailed and constructive review. 



APPENDIX 

A. ANALYTIC APPROXIMATIONS OF THE DISK-LIMITED GAS ACCRETION RATES 

In this section, we provide analytic approximations for the gas accretion rate in the regime where this rate is limited 
by the ability of the disk to transfer gas to the planet. In the calculations, we used piece-wise functions obtained by 
fitting the data from the 3D hydrodynamical calculations (see Section 3), for various values of the turbulent viscosity 
parameter, at, which quantifies the kinematic viscosity of the disk at the radial location of the planet, v t — a t H p Q p . 
We recall that the hydrodynamical calculations used an aspect ratio H p /a p = 0.05, which is a reasonable value in 
evolved disks between 5 and 10 AU (e.g., D'Angelo & Marzari 2012, and references therein). 

We fitted (log M p , log M p ) data using multiple second-odcr polynomials, which were then smoothly joined in overlap- 
ping regions. Since this procedure is somewhat cumbersome, here we provide simpler analytic approximations derived 
from data in the range of M p /M+ from 10~ 4 to 10~ 2 . In the calculations, as explained in Section 3, disk-limited 
accretion sets in when M p > 0.25 Mj up . 

Let us introduce the four functions 



where q — 
a t = 10- 2 
used: 



h(q)=ao ■+ 
.f2(q) = b + 
/3(g) = c + 

fi(q)=do -+ 

Mp/M* and all logarithms are in base 10 



01 logg + a 2 (logg) 2 
61 logq + b 2 (log qf 
ci logg + c 2 (logg) 2 
d\ logg + d 2 (log(7) 2 , 
The coefficients a, , hi, Ci , and di 



(Al) 
(A2) 
(A3) 
(A4) 

are given in Table Al. For 



the following analytic approximation for the disk-limited gas accretion rate, in units of a p T, p fl pi may be 

/ h(q) if q < 0.001197 



h{q) 



if 

otherwise. 



logM p = 

For at — 4 x 10~ 3 , the following analytic approximation may be applied 

logM p = min [f 3 (q),fa(q)]. 



(A5) 



(A6) 



The fitting functions are shown in the upper panel of Figure Al, along with the data obtained from the 3D hydrody- 
namical calculations. 

In the range of the turbulent parameter at that we explored (4 x 10~ 4 < at < 0.02), the maximum of M p occurs at a 
ratio M p /M+ similar (within a factor of « 2) to the square root of the right-hand side of Equation (7), i.e., before gas 
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Table Al 

Coefficients in Equations (A1)-(A4) 





i = 


i = 1 


i = 2 


Gk 


-6.8345179 


-2.5152600 


-0.354296 


bi 


-12.471200 


-6.3732500 


-1.014440 


Ci 


-11.429400 


-5.0171300 


-0.697484 


di 


-17.819292 


-8.5200885 


-1.172181 




0.0001 0.0010 0.0100 0.1000 

M p /M, 




E re( (g/cm 2 ) 

Figure Al. Upper panel: disk-limited gas accretion rates versus the planet-to-star mass ratio (see also Figure 2). The dashed line is 
Equation (A5), the solid line is Equation (A6), and the symbols represent the simulations' data for the disk turbulent parameter at = 10 -2 
(open circles) and 4 X 10 -3 (filled pentagons). Lower panel: final mass of a planet accreting at a disk-limited gas accretion rate in various 
situations: at = 4 X 10 — 3 (diamonds), at = 10 -2 (circles), a p = 5.2 AU (thin lines), a p = 9.5 AU (thick lines), M* = 1 Mq (open symbols), 
and M* = 2 Mq (filled symbols). See text for further details. 

begins to be depleted significantly because of the formation of the density gap. The maximum of M p can be compared 
to the (steady state) accretion rate through the disk in absence of the planet, 37ri/ t E p (Lynden-Bell & Pringle 1974), 
at the radial location of the planet. In units of a^Spfip, this accretion rate can be written as 3Trat(H p /a p ) 2 , giving 
ps 10~ 4 and 2.4 x 10~ 4 for a t = 0.004 and 0.01, respectively. As can be seen in Figure Al, these disk accretion rates are 
smaller than the maximum of M p . However, it should be noted that the tidal perturbation of the planet can modify 
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the accretion rate through the disk (Lubow & D'Angelo 2006). 

Equations (A5) and (A6) can be integrated to find the final (asymptotic) mass of a planet, M nna i, that accretes gas 
at a disk-limited gas accretion rate. We solved numerically the differential equation 

Alp = F(M P , M*, a p , E p , a t ) (A7) 

for M p , using an adaptive Adams-Bashforth-Moulton method of variable order with adaptive step-size and error 
control (available through the SLATEC Common Mathematical Library). Notice from Equation (A7) that, although 
the dependence of M p on M*, a p , and E p is trivial, the dependence of M p (t) on those three quantities is not! 

During the integration of Equation (A7), we assumed that M*, a p , and at are constants. In oder to mimic the 
viscous evolution of the (unperturbed) gas surface density at the radial position of the planet, £„, we applied the 
solution of Lynden-Bcll & Pringlc (1974) for a disk with no central couple. Using the same notations and indicating 
with i?i the initial standard deviation of the (gaussian) surface density distribution and with Mi the initial disk mass, 
Lynden-Bell & Pringle (1974) found that v t /R\ = (2/3)M*/Mi (where M* is the initial accretion rate on the star). 
Introducing the non-dimensional 'viscous' time^ t v j s = d>{v t / R\)t + 1, which can be written as t v i s = 4(M*/Mi)f + 1, 
the surface density evolution can be approximated by 

E p =S ref ^ /4 , (A8) 

where E ro f is a parameter and which represents the behavior of the Lynden-Bell & Pringle solution for i V j S 3> 1. We 
assumed that M*/Mi is ps 10 _6 yr _1 . According to the equation above, the surface density ratio E p /S re f decreases 
by more than two orders of magnitude over 10 Myr. 

We integrated Equation (A7) for the values of M*, a p , and a t used in the calculations, applying Equation (A8), and 
determined Mg na i as a function of S re f. The results are shown in the lower panel of Figure Al (see figure caption for a 

description of the different curves). The final mass is reached within about 5.5 Myr, when typically M p /M p ~ 100 Myr. 
The effect of disk viscosity is evident in this figure. In fact, around a solar mass star, the mass threshold for deuterium 
burning can only be achieved for at > 10 -2 . Among the varied parameters, at produces the largest differences in 
Affinal, whereas a p produces the smallest. Notice that the values of M nna i shown in the lower panel of Figure Al 
should not necessarily agree with those in the D-burning calculations because of the different assumptions made for 
the nebula evolution. In particular, S p in those calculations was taken as a constant. 



t Notice that the power of R± , in the definition of Lyndcn-Boll & ( . .. . , . u ■ , .. . n 

_ . , /-,„_,> , \ . , nil i i ■ in ■ r> in, refers to the viscous time t v ; s = 1, when the physical time t = 0. 

Pringle (1974), should be -2. Also, the subscript '1' in Ri and Mi vis , v j 
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